; docformat = 'rst'
;
;+
;
; :Purpose:
;   Read MOZART output from GAUGE folders
;
; :Inputs:
;   species
;   run_name
;   variable_name
;   first_year
;
; :Keywords:
;
; :Outputs:
;
; :Example::
;   mozart_output=gauge_mozart_read_history('CO2', 'CO2', 'CO2_VMR_avrg', 2004)
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Apr 2, 2014
;
;-
function gauge_mozart_read_history, species, run_name, variable_name, first_year, last_year=last_year

  year=first_year  
  date=!null
  fnames=!null
  elementsPerFile=!null
  
  if ~keyword_set(last_year) then begin
    last_year=2016
  endif
  
  repeat begin
    fnames_year=file_search(strcompress('/data/shared/GAUGE/MOZART/' + species + '/' + $
      string(year) + '/output/' + run_name + '.mz4.h0.*', /remove_all))
    if fnames_year ne '' then begin
      foreach fname, fnames_year do begin
        date_file=mr_ncget('date', filename=fname)
        date=[date, date_file]
        elementsPerFile=[elementsPerFile, n_elements(date_file)]
      endforeach
      fnames=[fnames, fnames_year]
    endif
    year++
  endrep until (fnames_year eq '') or (year ge last_year)


  lon=mr_ncget('lon', filename=fnames[0])
  lat=mr_ncget('lat', filename=fnames[0])
  lev=mr_ncget('lev', filename=fnames[0])
  hyam=mr_ncget('hyam', filename=fnames[0])
  hybm=mr_ncget('hybm', filename=fnames[0])

  case 1 of
    strmatch(variable_name, '*VMR*'): var=fltarr(n_elements(lon), n_elements(lat), n_elements(lev), n_elements(date))
    strmatch(variable_name, '*SRF_EMIS*'): var=fltarr(n_elements(lon), n_elements(lat), n_elements(date))
    strmatch(variable_name, '*LIFETIME*'): var=fltarr(n_elements(lon), n_elements(lat), n_elements(date))
  endcase
  PS=fltarr(n_elements(lon), n_elements(lat), n_elements(date))
  
  iFirst=0L
  foreach fname, fnames, ti do begin
    iLast=total(elementsPerFile[0:ti])-1
    case (size(var))[0] of
      4: var[*, *, *, iFirst:iLast]=$
        mr_ncget(variable_name, filename=fname)
      3: var[*, *, iFirst:iLast]=$
        mr_ncget(variable_name, filename=fname)
    endcase
    PS[*, *, iFirst:iLast]=$
      mr_ncget('PS', filename=fname)
    iFirst=iLast+1
  endforeach
  
  return, {var:var, ps:ps, lon:lon, lat:lat, lev:lev, hyam:hyam, hybm:hybm, date:date}

end
